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The excess free energy of a molecular cluster is a key quantity in models of the nucleation of 
droplets from a metastable vapour phase; it is often viewed as the free energy arising from the 
presence of an interface between the two phases. We show how this quantity can be extracted from 
simulations of the mechanical disassembly of a cluster using guide particles in molecular dynamics. 

We disassemble clusters ranging in size from 5 to 27 argon-like Lennard-Jones atoms, thermalised at 
60 K, and obtain excess free energies, by means of the Jarzynski equality, that are consistent with 
previous studies. We only simulate the cluster of interest, in contrast to approaches that require a 
series of comparisons to be made between clusters differing in size by one molecule. We discuss the 
advantages and disadvantages of the scheme and how it might be applied to more complex systems. 

PACS numbers: 82.60.Nh, 64.60.Q-, 36.40.-c 


I. INTRODUCTION 

The formation of droplets from a metastable vapour 
phase is a commonplace event in nature, but so far it 
has resisted quantitative analysis, despite repeated at¬ 
tention M- The phenomenon plays a role in atmo¬ 
spheric aerosol and cloud formation [5, (i|. as well as in 
industrial processes 01 - Theoretical analysis often be¬ 
gins with the Becker-Doring equations (9j that describe 
changes in the populations of clusters of i molecules 
brought about by the processes of gain and loss of single 
molecules, or monomers. They take the form 

drii/dt = Pi-iUi-i + a i+ 1 n i+1 - (fa + at) nt , ( 1 ) 

where fa and at are growth and evaporation rates, re¬ 
spectively. The rate of nucleation J of droplets from a 
metastable vapour phase may then be expressed as 0 

J = nifa*Zexp [- (4>(i*) - </>(l)) /kT\ , (2) 

where k is the Boltzmann constant, T is the tempera¬ 
ture, i* is the size of the critical cluster, defined to have 
equal probabilities, per unit time, of molecular gain or 
loss, and Z is the Zeldovich factor that accounts for the 
nonequilibrium nature of the kinetics a We shall re¬ 
fer to <j>(i) as the thermodynamic work of formation of 
a cluster of i particles (or i-cluster) starting from the 
metastable vapour phase. A range of nomenclature is 
used for this quantity in the nucleation theory literature: 
the work of formation was denoted by e(i) in [lGj, and 
elsewhere the same, or a very similar quantity has been 
labelled as A F, AG or AW, for example. 

We note that 4> in the nucleation rate expression has 
both a kinetic and a thermodynamic interpretation [ 0 . 
The quantity <j>(i) — <f>( 1) can be expressed in terms of 
ratios of cluster growth and evaporation rates: 

4>{i) ~ 0(1) = -fcr^h 1 ^-, (3) 

a i 



Figure 1. Typical work of formation of a cluster of i particles 
with a maximum at the critical cluster size i*. 

but <p is also related to the grand potential H s (*) = F(i) — 
ijjL s of an Acluster at the chemical potential fx s of the 
saturated vapour 0 : 

<f>(i) = Ul s (i) — ikTlnS, (4) 

where F(i) is the Helmholtz free energy of the clus¬ 
ter, S = Pv/Pvs is the vapour supersaturation, and p v 
and p vs are the vapour pressure and saturated vapour 
pressure, respectively. The role of the grand poten¬ 
tial in this context is to specify the equilibrium popu¬ 
lation of clusters of size i in a saturated vapour, namely 
nf = exp (—Sl s (i)/kT). The nucleation model is com¬ 
pleted by representing the population of monomers as 
n\ = Sp vs V/kT, where V is the system volume, by as¬ 
suming that the vapour pressure is dominated by the 
ideal partial pressure of single molecules. 

In classical nucleation theory (CNT), clusters are 
viewed as scaled down versions of macroscopic droplets. 
According to this approach, the difference <p(i) — (f(l) is 
replaced by <t>(i) alone with 

<t>(i) ~ <t>d(i) = ^A(i) - ikTlnS, (5) 

where 7 is the surface tension of a planar interface be¬ 
tween vapour and condensate, and A(i) is the surface 
area of a cluster represented as a sphere with a density 
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equal to that of the bulk condensed phase. The work of 
formation is a combination of a free energy cost of form¬ 
ing the interface, and a free energy return proportional to 
the number of molecules in the cluster (or proportional 
to its volume since the condensed phase density is taken 
to be a constant). The neglected (f)(1) term might be rep¬ 
resented by 7^4(1) — kT In S, which leads to the internally 
consistent classical theory [13|. 

The cluster size dependence of the CNT work of forma¬ 
tion is illustrated in Figure [1] It represents a thermody¬ 
namic barrier, with a maximum at the critical size, that 
limits the natural tendency for small molecular clusters 
to grow into large droplets when exposed to a supersat¬ 
urated vapour. CNT has been modified in several ways, 
for example by introducing a size-dependent surface ten¬ 
sion 0 or by introducing compatibility with nonideal 
vapour properties 0, 0. 

More fundamentally, the ratio of kinetic coefficients 
fij-i/aj might be evaluated using an underlying micro¬ 
scopic model for all clusters up to the critical size and be¬ 
yond 0, and the work of formation determined through 
Eq. ®. It may be shown that 

Pj-i/ a j = s ex P - ^s(j ~ 1)]/A;T], (6) 

which shifts attention to the free energy difference F(j) — 
F(j — 1) associated with the addition of a molecule to a 
(j — l)-cluster. Computing these differences is the basis 
of an approach has been used extensively in calculations 
of cluster free energies fl7l - l23j . But nucleation is actu¬ 
ally controlled by the properties of clusters near the crit¬ 
ical size, and one drawback of computing the differences 
F(j) — F(j — 1) is that the predicted nucleation rate could 
be susceptible to the accumulation of errors in evaluating 
such a sequence. 

In this paper, we describe a computational method 
for directly obtaining the cluster free energy without the 
need to perform calculations for a sequence of smaller 
clusters. We consider the following representation of the 
work of formation of a cluster minus that of a monomer: 

<t>(i) - 0(1) = F s (i) - (i - l)kTlnS. (7) 

We shall refer to F s (i) as the cluster excess free energy, 
though more accurately it is a difference between the ex¬ 
cess free energies of an Acluster and a monomer 0. It 
is ‘excess’ in that it represents the free energy required to 
carve a cluster out of a bulk condensed phase, or equiv¬ 
alently to assemble it out of saturated vapour. It may 
be associated with the thermodynamic cost of creating 
an interface, which is why in CNT it is modelled by a 
surface term, and why we have given it a suffix s. 

Our approach centres on disassembling a cluster into 
its component molecules using guided molecular dynam¬ 
ics in order to calculate the cluster excess free energy 
directly. The method employs the Jarzynski equality 
|24h26I | and we provide details in Section El including a 
comparison with the related method of thermodynamic 
integration. Tests of the method where we separate a 
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Figure 2. Guided disassembly process for an i-cluster. The 
real particles (circles) are initially weakly tethered to the 
guide particles (diamonds). The latter drift apart and the 
tethers gradually tighten leading to i independent, tethered 
particles upon completion of the process. 


dimer according to a variety of protocols are described 
in Appendix [A] The disassembly of argon-like Lennard- 
Jones clusters is presented in Section ITTT1 and we compare 
our results with those obtained from Monte Carlo studies 
by Barrett and Knight j2?i] and Merikanto etal. [28:, ,0. 
These studies gave consistent excess free energies, though 
they were not in agreement with experiments by Hand 
etal. [30]. We conclude with a discussion of the advan¬ 
tages and disadvantages of the approach compared with 
other treatments in Section Hvl 


II. GUIDED MOLECULAR DYNAMICS 
SIMULATIONS 


A. Fundamentals of the method 

We study the dynamical evolution of a cluster against 
a background of external manipulation. The cluster par¬ 
ticles are harmonically tethered to a set of artificial ‘guide 
particles’, which lie initially at the origin but after a pe¬ 
riod of system equilibration are programmed to move 
apart, driving cluster disassembly. The strength of the 
tether forces is initially quite weak, in order to disturb 
the properties of the cluster as little as possible. Later, 
the tethers can be strengthened in order to guide the 
separation process more firmly, and to prevent the atoms 
from interacting with each other once the final guide par¬ 
ticle positions have been reached. The mechanical work 
of the disassembly can then be related to the change in 
Helmholtz free energy. 

The masses of the guide particles are taken to be very 
much greater than those of the cluster particles. This 
essentially fixes the trajectories of the guide particles in 
the molecular dynamics, in accordance with the velocities 
assigned to each at the beginning of the disassembly pro¬ 
cess. By choosing guide particle velocities, simulation 
times and a time-dependent tethering force, a range of 
cluster disassembly protocols can be explored. A simple 
illustration of the process is shown in Figure [5] 

We shall consider clusters of argon-like atoms interact¬ 
ing through Lennard-Jones potentials, and so we shall re¬ 
fer to the cluster particles as atoms. We equilibrate this 
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system under the influence of the tethers for a suitable 
period, the duration of which will depend upon the clus¬ 
ter size and the desired temperature. A further molecular 
dynamics simulation is performed and from this trajec¬ 
tory we select initial configurations for cluster disassem¬ 
bly. In order that the configurations should represent a 
bound structure, we employ a Stillinger cluster condition 
[3l| in the selection, allowing a separation of no more 
than 1.5(TArAr between an atom and its nearest neigh¬ 
bour, where CTArAr is the usual Lennard-Jones range pa¬ 
rameter. Such a Stillinger condition has been used in 
previous Monte Carlo approaches. The cluster definition 
is an important ingredient of a modelling strategy [2j, 
and deserves careful consideration, but here we shall use 
this simple criterion for convenience. 

The simulations were performed using the DL_POLY 
[ 3 ^ molecular dynamics package, with modifications to 
the source code to implement the time-dependent har¬ 
monic tether potentials. We include a physical heat bath 
of helium-like Lennard-Jones atoms thermalised using a 
Nose-Hoover thermostat [33|. We could instead have im¬ 
plemented a thermostat that acts on the cluster itself, but 
chose not to in order to achieve as natural a thermalisa- 
tion as possible during the nonequilibrium processing. 


B. Work performed on a system 


Given an external control parameter A in a Hamilto¬ 
nian H( A), the work W done on a system due to the 
evolution of A over a finite time period may be written 


W = 


dX dH{ A) 
dt dX 


dt. 


( 8 ) 


For example, consider the Hamiltonian Hi of a single 
guided atom of mass to: 

Hi = ^+ l -K{t)[x(t)-Mt)]\ (9) 

where p is the momentum, n(t) is the time-dependent 
tethering force or spring constant, x(t) is the atomic po¬ 
sition and X(f) is the guide position. For a set of guided 
atoms, each controlled by a Hamiltonian H containing 
terms of the form given in Eq. Q supplemented by in¬ 
terparticle interactions, n(t) and X(f) play the role of A 
and the work W performed on the set is 


v j dXjjt) dH(K,{X k }) ^ 

J dt dn ^ J dt <9X 7 

0 J 0 

= \] 

0 j =i 

r 

- / «(*) E [Xj{t) - Xj(t)} ■ Vj(t)dt, 

0 J =1 

( 10 ) 


where r is the length of the molecular dynamics simula¬ 
tion, and V, ( t ) is the velocity of the guide particle asso¬ 
ciated with the j th atom, defined as = dX.j(t)/dt. 

The first term in Eq. uni arises from the time depen¬ 
dence of the spring constant, and the second term is sim¬ 
ply the conventional force times distance expression. It 
should be noted that all tethers within the system are 
characterised by the same spring constant, although more 
elaborate protocols could be imagined. 


C. The Jarzynski equality 

If we were able to perform an extremely slow, qua¬ 
sistatic process, then the mechanical work done would 
be equal to the difference in Helmholtz free energy be¬ 
tween the initial and final equilibrium states. However, 
quasistatic processes are unfeasible in finite time molecu¬ 
lar dynamics simulations and according to the second law 
0, the average of the work done (as a result of a time- 
dependent change in the Hamiltonian of the system), per¬ 
formed over many realisations of a nonquasistatic process 
(indicated by angled brackets), will always be an overes¬ 
timate of the free energy change, ( W ) > AF, allowing us 
only to infer an upper limit to AF. 

However, the Jarzynski equality 0.10 

(exp (— W/kT )) = exp (—A F/kT) (11) 

allows us to do better. For this identity to hold, the 
system must begin in thermal equilibrium, but need not 
remain so as the Hamiltonian changes during the simula¬ 
tion. Exploiting the work done in a nonequilibrium pro¬ 
cess is a powerful strategy for calculating cluster surface 
free energies and numerous computational studies (jaMjI) 
as well as experiments [j2l (d~7j | have achieved this with the 
help of the Jarzynski equality. Systems studied include 
argon-like Lennard-Jones fluids, ion-charging in water, 
ideal gases confined to a piston, and one-dimensional 
polymer chains. Nevertheless, there are distinct aspects 
of this strategy for analysing the controlled disassembly 
of a cluster that need to be explored. 

The Jarzynski equality ought to recover the free en¬ 
ergy difference regardless of the nature of the evolution 
between initial and final Hamiltonians, but computed re¬ 
sults might still depend upon the rate of the process as 
a consequence of a limited sampling of system trajecto¬ 
ries in finite simulations [42;]. We might expect ‘slow’ 
processes that gently pull a cluster apart to generate a 
narrower distribution of work compared with ‘fast’ pro¬ 
cesses that are violent and highly dissipative. A balance 
must therefore be struck between the poorer convergence 
of fast simulations and the demand for computational re¬ 
sources required for slow simulations. 

Furthermore, a consequence of the exponential averag¬ 
ing in the Jarzynski equality is that occasional values of 
work that are well below the average, arising from un¬ 
usual trajectories, can sometimes distort the extracted 
free energy change. This is a consequence of insufficient 
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sampling of the system trajectories and so we need to 
give careful attention to the statistical errors. 

We have explored the outcomes of various guiding pro¬ 
tocols, and the robustness of the Jarzynski equality in the 
face of limited statistics, in a test case of the separation 
of a dimer, for which the free energy change is easily cal¬ 
culable. These studies are described in Appendix [X] We 
have used similar protocols to study the disassembly of 
larger clusters, which is described in Section Hill 

D. Comparison with thermodynamic integration 

The method bears some similarity to thermodynamic 
integration, where the strength of the interparticle inter¬ 
actions is evolved over a sequence of equilibrium calcu¬ 
lations in order to compare the system in question with 
another that has a known free energy [48H5ll| . The basic 
relationship A F = f (dH(X)/d\)d\ is analogous to Eq. 
©. The reference system for clusters might, for example, 
be a set of noninteracting particles held together through 
the retention of the constraining cluster definition. Or in¬ 
deed the cluster definition could be changed progressively 
along with the interactions in order to reach a more con¬ 
venient final state, perhaps noninteracting particles in¬ 
side a sphere. 

However, there are some important differences. In our 
approach it is the tether potentials that change with time, 
not the interparticle interactions, and our reference sys¬ 
tem is a set of independent harmonic oscillators, not an 
ideal gas. Furthermore, we conduct the evolution by 
nonequilibrium molecular dynamics rather than by mov¬ 
ing through a sequence of equilibrium ensembles, and we 
only need to impose a cluster definition when selecting 
the initial configurations, not throughout the evolution. 
An abrupt removal of the cluster definition constraint is 
acceptable in a nonequilibrium evolution, when the re¬ 
sults are processed using the Jarzynski equation, but it 
would not be appropriate during a sequence of equilib¬ 
rium calculations. 


III. ARGON CLUSTER DISASSEMBLY 
A. Preliminaries 

We have investigated the disassembly of clusters con¬ 
sisting of 5, 10, 15, 20 and 27 argon-like atoms in or¬ 
der to obtain their excess free energies. Scaling up the 
guided molecular dynamics simulations from the test case 
of dimer separation is fairly straightforward. We perform 
simulations in a cubic cell with edge lengths of 100 A, so 
that the initial clusters and the final disassembled con¬ 
figurations may be easily accommodated. We employ 
Lennard-Jones interaction potentials for each species (see 
Table ED and the helium temperature is set at 60 K in 
order to facilitate a comparison with the Monte Carlo 
studies by Barrett and Knight and Merikanto etal. 



Figure 3. The difference in tether energy across a cluster con¬ 
figuration is given in terms of the maximum and minimum 
separations between an atom and its guide particle. The cir¬ 
cles depict the argon atoms, while the diamond represents the 
position of all of the guides at the origin of the cell. 

2^, J22|, as well as the experimental studies of Hand et al. 

However, converting the free energy change associated 
with disassembly into an excess free energy requires some 
careful consideration of the statistical mechanics of teth¬ 
ered and free molecular clusters. We require the excess 
free energy of a cluster that is free to move anywhere 
inside a system volume, but our initial state is a clus¬ 
ter tethered to guide particles at the origin. The free 
energy change that emerges from our calculations will 
correspond to the disassembly of a cluster whose centre 
of mass explores a region around the origin, and further¬ 
more, one that possesses energy due to the tethers in 
addition to that of the physical interactions between the 
atoms. These matters are discussed in detail in Appendix 
B. 

The energetic perturbation of the cluster by the tethers 
can be reduced by choosing a small force constant. We 
take the view that the mean variation in tethering energy 
of an atom, as it explores different regions of the cluster 
during the equilibrated trajectory, should not exceed the 
thermal energy kT, or 

7} K i (^max ~ ^min) < T, (12) 

where x max and £ m i n are, respectively, the maximum and 
minimum separations between an atom and its guide par¬ 
ticle in a configuration (see Figure [31). This criterion may 
also be expressed as £ = Ki(^ ax — x^ in )/(2 kT) < 1. 

From the equilibrated molecular dynamics trajectory, 
we select, for disassembly, a set of ‘valid’ cluster config¬ 
urations that satisfy the Stillinger cluster definition [31l |. 
but this can be quite difficult for the smaller clusters at 60 
K. Tethering the atoms keeps them closer together and 
more likely to form valid configurations. We therefore 
choose a tethering strength that satisfies the condition 
on £, but also helps to produce sufficient valid cluster 
configurations. The initial value of the tethering force 
constant was taken to be n t = 0.01 kJ mol _1 A , which 
gives £ ~ 0.6 — 0.9 for the five sizes of argon cluster 
studied. Table |T] shows the duration of the equilibrated 
cluster trajectory, the number of valid cluster configura¬ 
tions identified from candidates selected at intervals of 
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i 

Duration/ns 

Valid configurations 

i 

5 

1000 

152 

0.604 

10 

250 

411 

0.745 

15 

225 

1070 

0.799 

20 

150 

905 

0.847 

27 

150 

1020 

0.922 


Table I. The duration of equilibrated cluster trajectories at 
60 K, as well as the number of valid cluster configurations se¬ 
lected at each size. The ratio £ characterises the perturbation 
to the cluster energy due to the tether potentials. 


100 ps from the equilibrated molecular dynamics trajec¬ 
tory, and the ratio £ characterising the suitability of the 
tethering force constant. 

Having obtained initial cluster configurations for the 
five sizes of cluster, the next stage is to disassemble them 
by a combination of guide particle motion and tether 
tightening. A range of separation times t sep is explored, 
with the larger and more stable clusters expected to re¬ 
quire longer disassembly processes in order to provide 
accurate estimates of the free energy change. As in the 
dimer calculations described in Appendix A, we use a 
tethering strength that strengthens in time according to 
Eqs. (IA7I) . with a final value of k/ = 0.05 kJ mol~ 1 A 2 . 

The terminal positions for the guide particles are cho¬ 
sen from a 3 x 3 x 3 grid with spacing of 33.33 A. The 
largest cluster considered contains 27 argon atoms so af¬ 
ter the process of disassembly, the tethered atoms move 
around each point on this grid. For smaller systems, the 
same grid of final guide positions is adopted, but employ¬ 
ing only as many points as are necessary for the cluster in 
question. With initial guide positions at the origin and 
final positions defined in this way, it is straightforward to 
calculate the necessary drift velocities of the guide par¬ 
ticles for a given separation time. Applying the Jarzyn- 
ski procedure to the distribution of performed work then 
gives us the estimated free energy change AF associated 
with the disassembly of a cluster. 

However, as mentioned previously, this free energy dif¬ 
ference will only correspond to the disassembly of a teth¬ 
ered i-cluster, rather than of a freely translating, undis¬ 
torted cluster. Furthermore, by necessity we obtain free 
energies of systems of distinguishable atoms in molecu¬ 
lar dynamics, and we need to make an indistinguisha- 
bility correction. An analysis of the thermodynamics is 
required in order to extract the excess free energy of an 
f-cluster from the free energy of disassembly, and the de¬ 
tails are given in Appendix B. It turns out that we can 
write F s (i) = ELi /«*(*) with 


-A F 

(13) 

—ikT In (jO„sUho) 

(14) 

kT In (p vs v c ) 

(15) 

3 ini ( 3ivi \ 2 ^ 3 

(16) 

To” vItT ) 

kT In i ! 

(17) 


In the first term the free energy of disassembly A F ap¬ 
pears with a negative sign because it refers to the pro¬ 
cess of taking a cluster apart while F s is the free energy 
of interface formation. The / 2 term arises from relat¬ 
ing the final state in the disassembly process, namely 
the separated harmonically bound particles, to the ap¬ 
propriate reference state of a saturated vapour. It rep¬ 
resents the difference in free energy between the teth¬ 
ered particles, each effectively confined to a volume 
uho = (2-7T kT/ Kf) 3 ^ 2 , and particles in the saturated 
vapour phase with density p vs and volume per particle 
1/pvs■ The term is the entropy penalty associated 
with the initial tethering: the centre of mass of the cluster 
is effectively confined to a volume v c = (ini/(2irkT)) 3 ^ 2 
and needs to be referred to a situation where it is allowed, 
like a particle in saturated vapour, to explore a volume 
1 /p vs . The term is an approximate expression for 
the perturbation in the cluster energy due to the initial 
presence of the tethers, where V[ = 1/pi is the volume 
per particle in the condensed phase. Finally, f% converts 
calculations derived from molecular dynamics with dis¬ 
tinguishable particles into results relevant to a system of 
indistinguishable particles. 


B. Results and discussion 

A typical example of the work W ( t ) performed over 
a disassembly trajectory of duration 20 ns for a 27-atom 
cluster is shown in Figure [4] The gradual rise in the work 
performed prior to about 5 ns represents an accumulation 
of tethering energy as the guide particles move away from 
their initial positions at the origin. After this time, atoms 
begin to leave the cluster, and less work is needed to move 
the corresponding guides. After about 7 ns, the work rate 
reduces significantly as the cluster disintegrates and the 
guide particles move towards their final positions. 

Visual representations of the disassembly process (see 
Figure 0 provide further insight into the manner in 
which the clusters are pulled apart. The onset of cluster 
disassembly is signalled by the loss of one or two atoms 
from the cluster, perhaps only temporarily. The clus¬ 
ter soon after breaks into several smaller clusters, which 
eventually disintegrate into fragments or single atoms. It 
is rare to see a complete and sudden disintegration of a 
cluster, where all the constituent atoms disassemble to¬ 
gether within a short space of time. 

Figures [6] and [7] show distributions of the work per¬ 
formed in disassembling the 5-cluster and the 27-cluster, 
along with estimates of the free energy change, for sepa¬ 
ration times between 0.5 ns and 20 ns. As expected, the 
work distributions are broader for the processes that are 
most rapid (smallest t“ p ) and hence least quasistatic in 
nature. Conversely, the work distributions become nar¬ 
rower, and lead to free energy changes that presumably 
provide the most accurate estimates of the true free en¬ 
ergy change, as the rate of separation is reduced. 

The free energy change A F for the disassembly of each 
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Figure 4. A typical history of the work performed for one 
realisation of the disassembly of a 27-atom argon cluster with 
a separation time of 20 ns. 


i 

tsep(ns) 

(W) 

A F 

/.(0 

fs(i) 

/»(*) 

fHi) 

F s (i) 

5 

6 

13.08 

12.35 

38.94 

-7.79 

-0.41 

4.79 

23.18 

10 

8 

34.06 

30.80 

77.88 

-8.83 

-1.32 

15.10 

52.03 

15 

12 

62.75 

53.87 

116.81 

-9.44 

-2.59 

27.90 

78.82 

20 

16 

97.37 

84.07 

155.75 

-9.87 

-4.18 

42.33 

99.97 

27 

20 

154.41 

133.28 

210.27 

-10.32 

-6.90 

64.56 

124.33 


Table II. Results from the slowest set of disassembly simu¬ 
lations for each cluster size: the mean work ( W), the free 
energy of disassembly A F and the other contributions to the 
excess free energy F„(i), all in units of kT. 


size of cluster at the slowest rate studied is shown in 
Table [Tl] along with the other contributions to the ex¬ 
cess free energy F s . We refer to a molecular dynamics 
study by Baidakov et al. [53| to provide values of the sat¬ 
urated vapour density p vs and liquid density pi = 1/vi 
of the argon-like Lennard-Jones fluid at a temperature of 
60.31 K. 

Figure [5] shows our excess free energies F s (i) as a func¬ 
tion of cluster size i. Statistical errors propagated from 
uncertainties in the free energy change A F are similar 
to the size of the symbols. We also include correspond¬ 
ing results from the Monte Carlo studies by Barrett and 
Knight [23] and Merikanto etal. [28l [29]. Barrett and 
Knight employed a Lee-Barker-Abraham cluster defini¬ 
tion E3 while Merikanto etal. adopted a Stillinger clus¬ 
ter criterion similar to ours. The Barrett and Knight cal¬ 
culations are represented here by Ff K (i)/kT = — In (p — 
(i — 1 ) ln(p„ s cr| rAr ) with their fitting function lng^ = 
10.5 + 9.91(? — 1 ) — 16.36(i 2 / 3 — 1), and the Merikanto 
et al. values are derived from their Figure 1 in (Hj, which 
we interpret as a plot of F^{i)/kT — (i — l)lnS with 
S = 20. The results of these earlier studies are consis¬ 
tent with one another, as well as with the excess free 
energy suggested by the internally consistent classical 
theory (ICCT) F^ CGT (i) = 7 (36nvf) 1/3 (i 2/3 - l) 



Figure 5. Illustration of the disassembly of a 27-atom ar¬ 
gon cluster, with green spheres representing the argon atoms 
and lighter spheres the guide particles (helium atoms are not 
shown). In frame 1, all the guides lie at the origin of the 
cell. By frame 2, the guides have drifted far enough apart for 
a single argon atom to escape temporarily from the cluster 
before rejoining it in frame 3. In frame 4, several atoms have 
escaped, but remain in close proximity to the reduced cluster. 
A threshold is reached in frame 5, where many argon atoms 
break free to leave a fragment of about five atoms that also 
soon disintegrates as shown in frame 6. Shortly after, all of 
the atoms fall into motion about their partner guide particles 
which continue along steady paths away from one another 
(frames 7 and 8). The reader is encouraged to view movies of 
the disassembly provided in the Supplemental Material [52 ], 
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Figure 6. The distribution of work W (top) for sets of disas¬ 
sembly trajectories for the 5-atom argon cluster, for a range 
of separation times. The lower plot shows the mean of the 
work l/W) and the corresponding free energy differences A F 
calculated via the Jarzynski equality for each t Bep . 


structed such that i 7 'j CCT (l) = 0 , where 7 is the surface 
tension of the planar liquid-vapour interface, again taken 
from Baidakov etal. [53j|. It is clear from Figure [5] that 
the calculations presented in this study are consistent 
with the previous Monte Carlo results. This is satisfac¬ 
tory support for the disassembly approach that we have 
developed. We note that all three are reasonably well 
represented by the ICCT model, which is somewhat sur¬ 
prising. 

Note that the construction of a traditional plot of the 
nucleation barrier such as Figure Q] would require us to 
subtract a term ikT In S from the excess free energies in 
Figure [ 8 ] Inserting a supersaturation of 30 would then 
yield a critical size of about 20 , for example. 


Figure 7. Plots similar to those shown in Figure HJ] but for the 
27-atom argon cluster. 



Figure 8. Excess free energies for argon-like Lennard-Jones 
clusters obtained from disassembly at 60 K are shown as 
squares and compared with values obtained in Monte Carlo 
studies by Barrett and Knight (27 | at 59.88 K (solid line) and 
Merikanto etal. [28i.[29] at 60.18K (triangles). Also shown is 
the prediction from internally consistent classical nucleation 
theory for a temperature of 60.31 K (dashed line). 
























IV. CONCLUSIONS 


We have developed a method of guided cluster dis¬ 
assembly in molecular dynamics, capable of extracting 
the excess free energy associated with the formation of a 
molecular cluster from the saturated vapour phase. This 
property is often regarded as a surface term and it plays 
a central role in kinetic and thermodynamic models of 
the process of droplet nucleation. 

After exploring some aspects of the method by separat¬ 
ing a dimer, the technique was applied to the controlled 
disassembly of Lennard-Jones argon clusters between 5 
and 27 atoms in size. The extracted free energy of disas¬ 
sembly has been related to the excess free energy of the 
cluster through an analysis of the statistical mechanics of 
free and tethered clusters. Our calculations for clusters 
of various sizes are consistent with previous studies by 
Barrett and Knight [271] and Merikanto et al. [28, 29], 
both of which require the evaluation of a sequence of free 
energy differences between monomer and dimer, dimer 
and trimer, etc. A Lennard-Jones microscopic model of 
argon, within the standard kinetic and thermodynamic 
framework of nucleation theory, cannot account for the 
experimental argon nucleation data of Hand et al. [3p|, 
but we do not speculate here about this disparity. 

The approach should be contrasted with methods of 
free energy estimation based on thermodynamic integra¬ 
tion. In those methods, the strength of the interparticle 
interactions is evolved over a sequence of equilibrium cal¬ 
culations. Our approach also involves the evolution of a 
Hamiltonian, but it is the tether potentials that change 
with time, not the interparticle interactions. Further¬ 
more, we evolve by nonequilibrium molecular dynamics 
rather than studying a sequence of equilibrium ensem¬ 
bles, and we are only required to apply a cluster defini¬ 
tion when selecting the initial configurations, not during 
the evolution. 

We believe that our process of mechanical disassem¬ 
bly offers an intuitive understanding of the meaning of 
the work of formation that plays such a central role in 
nucleation theory. We suggest that a direct evaluation 
of this quantity is preferable to an approach based on 
summing the free energy changes associated with the ad¬ 
dition of single molecules to a cluster, on the grounds that 
we avoid the possible compounding of statistical errors. 
The computational costs of our current study of argon 
clusters have been higher than those of more traditional 
methods such as grand canonical Monte Carlo [29], for 
the same level of accuracy, largely because of our explo¬ 
ration of different protocols and our use of an explicit 
helium thermostat, but these can be reduced with fur¬ 
ther development. A particularly powerful variant of the 
disassembly scheme is to separate a cluster into two sub¬ 
clusters under similar mechanical guidance, in order to 
relate the distribution of work performed to a free energy 
of ‘mitosis’, essentially a difference in excess free ener¬ 
gies between the initial cluster and the two final subclus¬ 
ters. Such comparisons would be unfeasible to perform 


in Monte Carlo. The calculations are not onerous and an 
evaluation of the excess free energy of clusters of up to 
128 water molecules is to be reported (54j. Furthermore, 
the explicit thermostat can be replaced by an implicit 
scheme. With such tools, and guided by the experience 
developed in the current investigation of argon, we in¬ 
tend to carry out studies of clusters of water, acids and 
organic molecules, species that are particularly relevant 
to the process of aerosol nucleation in the atmosphere. 
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Appendix A: Argon dimer separation 


We test the feasibility of the approach using two proto¬ 
cols of controlled dimer separation. First, the guide par¬ 
ticles are made to drift apart with the tether strengths 
held constant, and then we allow the tethers to tighten 
over the course of the process. We determine the maimer 
of dimer separation that leads to an accurate estimate of 
the free energy change. 

We start by evaluating the free energy of a tethered 
dimer of argon-like atoms analytically. Particles are 
distinguishable in molecular dynamics simulations since 
they carry labels, so we take this into account in the 
analysis. The initial Hamiltonian of the dimer system is 


^dimer 


2 2 
Pi ( Pi 

2 m 2 to 


( x 1 - Xi ) 2 


1 2 

-m{ aj2-X 2 ) +$(|®! -® 2 |) : 


(Al) 


where m is the argon mass, and <F(|£ei — a^l) is a pairwise 
interaction potential. When the guide particles both lie 
at the origin (Xi = X 2 = 0), the initial partition function 
is 


Z 


dimer 


= 7T / ex P ~ 


h e 


Pi 


■pI 


exp —m 


x\ 


x\ 


2 kT 


2mkT 
exp I - 


dp\dp 2 

■F ([ag - x 2 \) 

kT 


(A2) 
dxidx 2 , 


noting that there is no correction factor of one half since 
the atoms are distinguishable. Substituting r = X\ — x 2 
and R = xi + x 2 , the partition function Zf' mel becomes 


^th 


■ exp -Ki 


A t 6 h 


- / exp l -m 


1 IT f 4-7T kT\ 2 f 

^th ^ \ K i ) Jo 


2 kT 
R 2 + r 2 


exp - 


<F (|£Ci - x 2 \) 


AkT 


exp - 


kT 
<F (r) 


drdR 


r exp l — 


Ki.r 


kT 

2 - 1 - 4<F(7') 


r 2 drdR 


AkT 


dr , (A3) 
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where Ath = h/ (27rmfcT) 1 / 2 is the thermal de Broglie 
wavelength. We have imposed an upper limit r c on the 
separation between the two atoms, corresponding to a 
definition of what we mean by a dimer. 

For the final state in which the two argon atoms are 
tethered to respective guide particles that are far apart, 
the Hamiltonian is simply that in Eq. (EB without the 
interaction term, and with a final tether strength Kf. The 
corresponding final partition function is 


^dimer 


1 



Pi + v\ \ 

2 mkT J 


dp\dp 2 


(A4) 


X 


exp 



dx\dx2 


1 / 2irkT A 3 

A t 6 h V «/ ) 


The free energy change in separating a dimer of tethered 
atoms can therefore be expressed as 


A F = kT In (Zf mel /Zf mel ) 


= kT In 



(AS) 

K,;r 2 + 4$(r) A , 
- 4 W - ) dT ’ 


j k 

€jk / kJ mol 1 

o jk / A 

Ar Ar 

0.995581 

3.405 

He He 

0.084311 

2.600 

Ar He 

0.289721 

3.000 


Table III. Parameters for the Lennard-Jones potentials, where 
j and k are the atomic labels, tjk is the depth of the potential 
well, and <jjk is the range parameter |55j. 



Figure 9. Illustration of the dimer separation process. Both 
guide particles (diamonds) are initially at the origin, but one 
is made to drift towards a corner of the simulation cell. 


a. Guiding at constant tether strength 


which can be evaluated numerically. The parameter r c is 
the Stillinger radius used to identify a dimer configura¬ 
tion in the equilibrated molecular dynamics simulation, 
to which we now turn. 

We place two argon-like particles within a periodic cell 
with edge length 50 A, each tethered to guide particles 
through a harmonic interaction ^K(t)i' 2 , where r is the 
separation between the argon atom and its guide, and 
n{t) is the tethering force constant. The argon atoms 
are thermalised through interaction with a gas of 100 
helium-like atoms kept at constant temperature using a 
Nose-Hoover thermostat. Conventional masses of 39.85 
and 4.003 amu for the argon and helium-like particles are 
adopted, while the guide particles are assigned a vastly 
greater mass of 4 x 10 12 amu. Interaction potentials are 
specified by 


One of the guide particles drifts from the origin to a 
corner of the cubic simulation cell over a separation time 
f sep while the other remains stationary (see Figure EJ. 
We choose t sep to be 1, 2 or 4 ns and the velocity of 
the moving guide particle (labelled 1) is given by Vi = 
[Xi(t = tsep) — Xi(t = 0)] / t sep ■ 

For initial and final tethering force constants of 
0.05 kJ moP 1 A , the expected free energy change in 
separating the dimer is 5.716 kT according to Eq. Ell- 
Distributions of the work done for each rate of dimer 
separation are shown in Figure 1101 and the correspond¬ 
ing estimates of the free energy change obtained from the 
Jarzynski equality are compared with the expected value 
in the lower part of Figure fill A longer separation time 
leads to a better estimate of the free energy change since 
the process is then closer to being quasistatic. 


$ (pfc) 




(A6) 


with parameters shown in Table mu though it should be 
noted that only the repulsive part of the interaction be¬ 
tween argon and helium is employed in order to prevent 
any binding between the two. Simulations are performed 
at a temperature of 15 K such that dimers are long-lived 
and a sufficient number of configurations satisfying the 
separation criterion r < r c = 1.5<TArAr can be obtained 
from the equilibrated trajectory. With a constant tether¬ 
ing force constant of 0.05 kJ moP 1 A , we generate an 
equilibrated molecular dynamics trajectory of duration 
100 ns and choose 10 3 dimer configurations for use as 
starting points for the separation process. 


b. Guiding with tether tightening 

We now elaborate the process by tightening the tethers 
during guide drift according to 


n(t) = Ki for t <ti 

Kf - Ki 

- Ki + ^— 
= Kf for t > t s 


1 — cos 



for ti < t < t s 
(A7) 


where ti is the time at which the force constant begins 
to change, and t s is the time at which it reaches its final 
value. Once again starting with dimer configurations and 
an initial tethering force constant of 0.05 kJ moP 1 A 
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Figure 10. Distributions of the work done in the disassembly 
of a dimer for separation times t sep of 1, 2 and 4 ns. 


Figure 12. Distributions of the work of dimer disassembly 
where the atoms are guided apart and the tethers tightened 
for three different separation times. 
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Figure 11. Convergence of the Jarzynski-estimated free en¬ 
ergy change toward the expected value (dashed line) as the 
dimer separation rate is decreased, while keeping the tether¬ 
ing strength constant (lower set) and when the tethers are 
tightened (upper set). 


Appendix B: Analysis of cluster free energies 


1. Free and tethered clusters 


The canonical partition function Zp = exp (— Fp/kT ) 
for an untethered, or ‘free’ cluster of i indistinguishable 
particles governed by a Hamiltonian H composed of ki¬ 
netic energy terms and pairwise interactions is given by 


Z F = 


l\h 3i 


dxjdpj exp [-H ({a: fc }) /kT] ,(B1) 

i=i 


where Fp is the associated free energy. For a cluster 
tethered to the origin, the Hamiltonian will include an 
additional set of harmonic potentials, such that the par¬ 
tition function is 


at 15 K, three dimer separation times are investigated, 
during which the force constant rises by a factor of two. 
The times ti and t s are specified as 20% and 80% of the 
total separation time. The expected free energy change 
associated with dimer separation is 7.795 kT according 
to Eq. (IA5I) . It can be seen from the upper part of 
Figure [IT] that all three separation rates give acceptable 
estimates of the free energy change. Furthermore, the 
greater compatibility between the distributions of the 
work performed at different separation rates shown in 
Figure [T^I compared with those in the simulations with 
constant tether strength, suggests that a protocol where 
the tethers tighten while the guide particles drift apart 
is more effective. Intuitively, the separation is then con¬ 
ducted more firmly, and with less dissipation. 


Zp = exp (— Fp/kT ) 



x exp 




1=1 



/kT 


,(B2) 


where Fp is the free energy of the tethered cluster, and 
Ki is the initial tethering force constant. 

We insert a factor of unity in the form 1 = 

/ 5 l x j ~ x cj dx c into Eqs. (IBID and (IB2D . and 

transform to particle coordinates with respect to the clus¬ 
ter centre of mass x c , namely x/ = xj—x c . The partition 
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function for a free cluster becomes 



where V is the system volume and Zf is the partition 
function for a cluster whose centre of mass is fixed at 
the origin. It should be noted that since the Hamilto¬ 
nian contains pairwise interactions, it may be rewritten 
as H({x k }) = H({x' k }) after the change of variables. 

Similarly, the partition function for a tethered cluster 
can be rewritten as 


Zp = 


x exp 


Y[dx'jd Pj dx c 5 7 I>; (B4) 


l=i 


l=i 


i\h 3x 

- (^(K}) + E« 2 ) ! kT 


1=1 


The second term in the exponent of Eq. (IB4I) may be 
simplified using the constraint Xq=i x 'j = 0 an< 4 ^ follows 
that jy j= i Xj = Ej =i x f + ix % giving 


1 f /■ 

Zr= iW I Ii dx 'j d Pj dx c^P 


1=1 


--Hiixl/kT 


x exp 


- H ({x' k }) + /kT 


l=i 




1=1 


( o 7 rlT\ 2 i r * 

—-) ^3? I n dx 'j d Pj ex P \~ H (Kl) /kT\ 

1 ' ’ 1=1 


x exp 


-\^ii x V kT 


1=1 


4 tE-S 


1=i 


V IKi ) 


(B5) 


exp (—Fp/kT) = f dT exp [— (Hq + U) /kT], where T 
represents the configuration of a system, and dT is pro¬ 
portional to the phase space volume element Hjdx'jdpj. 
In the context of the tethered cluster described by Eq. 
(IB5I) . U represents the term Ylj=i x f> while H 0 is 
the untethered Hamiltonian H ({x' k }) modified by the 
delta function constraint. Fp is therefore the free energy 
of a tethered cluster with its centre of mass further con¬ 
strained to lie at the origin, and is equal to — kTlnZf. 
A similar relationship exists between Ff , the free energy 
of an untethered cluster with fixed centre of mass, and 

Zp- 

The free energies Fp and Fp may be related through 


exp (— Fp/kT ) = 


f dT exp (— Ho/kT) exp (—U/kT) 


f dT exp (—Ho/kT) 
x J dr exp (—Ho/kT) 

= (exp ( -U/kT)) 0 exp (-Ff/kT) ,(B6) 


where angle brackets represent an average in the statis¬ 
tical ensemble corresponding to Hq. For small ( U/kT ) 0 , 
we can write (exp (— U/kT)) 0 ~ exp (— (U) 0 /kT), and 
hence 


exp (-Ff/kT) ~ exp [(-F£ - (U) 0 ) / kT }, (B7) 


with (U )o given by 


_ / dTU ({x' k }) exp (-H 0 /kT) 
[ >0 fdr exp {-H 0 /kT) 


(B8) 


U({x' k }) is a sum of single-particle harmonic poten¬ 
tials of the form Uno( x k) = \ K i x k i so Eq. (IB8I) can be 
written as 


(U) o 


ELJdrt/HQ (x' k ) exp (—Hp/kT) 


f dr exp (—Ho/kT) 

. J dTUno (x' k ) exp (-Hp/kT) 
f dr exp (—Hp/kT) 


*(^ho)o-(B9) 


We next introduce the spatial density profile of a single 
particle (labelled k without loss of generality) in a cluster 
constrained to have its centre of mass at the origin but 
not tethered, namely 


. _ / dT exp (—Hp/kT) 5(x' k — y) 
P0[V> fdr exp (-Hp/kT) 


with f pp(y)dy = 1. We can write 


(BIO) 


where Zf is the partition function of a cluster constrained 
to have its centre of mass at the origin as well as having its 
constituent particles tethered to the origin by a harmonic 
potential. 

Next, we employ the Gibbs-Bogoliubov approach |56|, 
it] to compare the free energies Fp and Ff of sys¬ 
tems with Hamiltonians Hp and Hamiltonian Hp + U, 
defined by exp (—Fp/kT) = f dT exp [—Hp/kT] and 


(^ho)o = J Po(y)U ii o(y)dy 1 (Bll) 

which represents the average tethering energy of a parti¬ 
cle that is spatially distributed according to the density 
Po(y). The condition that the tether potential makes a 
relatively small contribution to the mean energy of the 
cluster is (Hho)o = / Po(y)y 2 dy <C kT, in which case 
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the approximations involved in the Gibbs-Bogoliubov ap¬ 
proach are acceptable and the initial tethering potential 
weak enough that the cluster is only slightly distorted in 
comparison with a free cluster. Thus we write 


Zp = exp (-Fp/kT) ~ exp [(-Fp - i (f7 H o) 0 ) /kT]. 

(B12) 


Eq. (IB5jl can then be written as 


Zp 


x exp 


n dx 'j d Pj ex p \~ h (md / fcr ] 

i=i 


/ 2nkT\ 2 1 

\ ini ) i\h 3x 

-i J Po (y) Kiy 2 dy/2kT S | j ^ s'- J , (B13) 


j =i 


such that the relationship between the partition function 
of a tethered cluster, and the partition function of a free 
cluster with a constrained centre of mass Zp, is 


Z T = Zj 


(2-KkT\ 2 
V *«» / 


exp 


-i J Po (y) my 2 dy/2kT 


Combining Eqs. drnli and dm then gives 
In Zp = In 


(B14) 


V \ IKi 


Zrf^kiy _ J My2dyA B15) 


— Fj 1 = — 


fcTln[pc(0)E] - ^ J p 0 {y)y 2 dy, (B16) 


where (iKi/2irkT) 3 ^ 2 has been replaced by a function 
p c ( 0), representing the probability density that the cen¬ 
tre of mass of the tethered cluster lies at the origin. This 
equivalence can be demonstrated by deriving the distri¬ 
bution of the cluster centre of mass, through considering 
a single particle with mass M = im and coordinates x c 
and p c residing in a potential iKiX 2 /2. The positional 
probability density at z is 


Pc(z) = 


f dx c dp c ex.p 


2 kT 


2 MkT 


) 6 ( X c - z) 


f dx c dp c exp 


= ( 


lKiX~ 

2 kT 
2 N 


2 MkT 


) 


Vssf) exp ( _ w) 


(B17) 


such that p c (0) = (iKi/(2nkT)) 3 ^ 2 . 

The purpose of the substitution is that the first term 
on the right hand side in Eq. dEU) may be interpreted as 
two competing contributions to the free energy difference 
Fp — Fp. We write 


— kT In \p c (0) V] = —T 



fclnE 


(B18) 

such that the first term corresponds to the removal of the 
entropic contribution to free energy associated with the 


freedom of motion of the cluster centre of mass within a 
constrained volume l/p c (0), brought about by the teth¬ 
ers, and the second term represents the addition of en¬ 
tropic free energy corresponding to the freedom of motion 
in volume V. Finally, the second term in Eq. (IB16I) is 
an estimate of the removal of tethering potential energy 
when relating a tethered to a free cluster. 


2. Excess free energy from the free energy of 
disassembly 

We now establish the relationship between the free 
energy of a free cluster to the cluster work of forma¬ 
tion defined as 4>(i) = ~ ikTlnS , where f = 

Fp(i) — ip s is the grand potential of a free cluster of 
i particles in an environment at chemical potential p s 
for which the bulk condensed and vapour phases coex¬ 
ist. The excess free energy (difference) of the cluster is 
therefore 


F s (?) = (f> (i) — 0(1) + (i — 1) kT In S 

= F F (i)-F(l)-(i-l)p s , (B19) 

having used Eq. (IB- 

Assuming the vapour is ideal, the coexistence chem¬ 
ical potential p s and the monomer Helmholtz free en¬ 
ergy F(l) are simply p s = kT\n(p vs A) and F(l) = 
—kT ln(V/A), respectively, where p vs is the particle den¬ 
sity in a saturated vapour and A = Aj? h with Ath = 
h/ (27rmfcT) 1 / 2 . The excess free energy F s (i) can now 
be expressed as 

F s {i) = Fp + kT In (V/A) - (i - 1) kT In (p vs A) 

= F t - kT In \p c (0) V] - -y- J p 0 (y) y 2 dy 
+kT In (V/A) - (i - 1) kT In (p vs A). (B20) 

Now we consider the free energy change associated 
with the process of cluster disassembly. The difference in 
free energy between separated constituent particles each 
tethered to a guide particle, and a tethered cluster, is 
SF = Ff — Fp , where Ff = —3ikT\n (kT/hujf ) is the 
free energy of i harmonic oscillators in three dimensions, 
where the angular frequency uif = ( Kf/m ) 7 of the os¬ 
cillators is related to the final value of the tethering force 
constant . 

It should be recognised, however, that the quantity 
SF is not the free energy difference extracted from the 
molecular dynamics simulations of cluster disassembly. 
Molecular dynamics simulations always involve distin¬ 
guishable particles, since they are assigned labels, and 
SF is a difference between the free energy of i indistin¬ 
guishable particles in a cluster, and i particles that are 
distinguishable through having been physically separated 
to regions around their final tether points. 

The free energy difference that is extracted in our pro¬ 
cedure is actually A F = Ff—Fp lst , where the superscript 
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in F^ lst reminds us that it is the free energy of a tethered 
cluster of distinguishable particles. But we can relate the 
partition function of such a cluster to the partition func¬ 
tion Zt for indistinguishable particles by the usual clas¬ 
sical procedure, namely Z^ lst = HZt , and since F^' st = 
—kT In Z^ lst = -kT hi Z T - kTlnil = F T — kT In i\ we 
have 


A F = F f -F t + kT In il = SF + kTlnil, (B21) 

such that Ft = Ff — SF = Ff — A F + kTlnil. Substi¬ 
tuting into Eq. (IB20I) then gives 

F s {i ) = — A F — ikT In (/ 9 „ s uho ) + kTlnil 

- kT In J p 0 (y)y 2 dy, (B22) 


where uro = (27 rkT/ kj) 3 ^ 2 is a volume scale associated 
with the confinement of particles within the final har¬ 
monic tether potentials. It should be noted that the ex¬ 
cess free energy F s does not depend upon the Planck 


constant h , nor on the system volume V, as is to be ex¬ 
pected. 

In order to complete our specification of F s (i) in terms 
of A F and material properties, we need to estimate 
the final term in Eq. (IB22I) . We write f po(y)y 2 dy = 
Jo°° Po(f)47rr 4 dr, where r is the distance from the cluster 
centre of mass, and recall that po(r) is the single-particle 
density profile in an untethered cluster with fixed centre 
of mass. As an approximation, we imagine the cluster to 
be spherical with a constant particle density, such that 
Po(r) ~ pi/i for 0 < r < r max , where pi is the parti¬ 
cle density in the condensed phase, and r max is the ra¬ 
dius of the cluster. Since the probability density po(r) 
is normalised, we have / 0 max (pi/*)47rr 2 dr = 1, such that 

r-max = (3i/4npi) 1/3 and so 


T max 



4?r Pi r 5 

5i max 


3 

5 



, (B23) 


where vi = l/pi is the volume per particle in the con¬ 
densed phase. Substituting this into Eq. (IB22II gives 
Eqs. (I13H17H in the main text. 
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